UNCLASSIFIED 


Defense  Technical  Information  Center 
Compilation  Part  Notice 

ADP013715 

TITLE:  Uniform  Powell-Sabin  Splines  for  the  Polygonal  Hole  Problem 
DISTRIBUTION:  Approved  for  public  release,  distribution  unlimited 


This  paper  is  part  of  the  following  report: 

TITLE:  Algorithms  For  Approximation  IV.  Proceedings  of  the  2001 
International  Symposium 

To  order  the  complete  compilation  report,  use:  ADA412833 

The  component  part  is  provided  here  to  allow  users  access  to  individually  authored  sections 
of  proceedings,  annals,  symposia,  etc.  However,  the  component  should  be  considered  within 
the  context  of  the  overall  compilation  report  and  not  as  a stand-alone  technical  report. 

The  following  component  part  numbers  comprise  the  compilation  report: 

ADP013708  thru  ADPO 13761 


UNCLASSIFIED 


Uniform  Powell-Sabin  splines  for  the  polygonal 

hole  problem 

Joris  Windmolders  and  Paul  Dierckx 

Department  of  Computer  Sciences,  Kath.  University  Leuven,  Belgium. 
Joris . WindmoldersQcs . kuleuven .ac.be,  Paul . Dier ckxOcs . kuleuven .ac.be 


Abstract 

An  algorithm  is  described  for  smoothly  filling  in  a polygonal  hole  in  a surface,  with  a 
parametric  uniform  Powell-Sabin  spline  surface  patch.  It  uses  interpolation  and  subdi- 
vision techniques  for  iteratively  determining  an  approximating  solution.  No  assumptions 
are  made  about  the  surrounding  surface.  The  user  has  to  provide  routines  for  calculating 
the  curve  points  and  the  unit  surface  normal  along  the  edge,  as  well  as  the  unit  tangent 
vector  of  the  edge  curves,  parametrized  on  the  unit  interval. 


1 Introduction 

A classical  problem  in  CAGD  is  to  fill  in  a hole,  bounded  by  a set  of  surfaces.  This 
problem  has  already  been  addressed  in  the  literature  (e.g.  [1,  2,  4]).  In  most  cases, 
assumptions  are  made  on  the  bounding  surfaces.  In  this  paper,  we  present  an  algorithm 
for  filling  in  a 3,  4,  5 or  6-sided  hole  that  makes  no  assumptions  on  the  surrounding 
surfaces,  and  therefore  it  is  generally  applicable.  On  the  other  hand,  the  filling  patch 
will  meet  the  given  boundary  curves  approximately.  The  input  of  our  algorithm  (see 
Figure  1)  consists  of  the  boundary  curves  p which  join  at  their  endpoints.  Furthermore, 
the  user  should  provide  the  unit  tangent  vector  7 to  the  boundary  curves  at  any  point, 
and  the  unit  normal  vector  n to  the  surrounding  surface  at  any  curve  point  except  the 
endpoints,  where  the  tangent  vectors  of  the  joining  curves  are  needed  only  (see  Figure 
1 again).  For  other  (interior)  curve  points,  our  algorithm  will  calculate  a unit  vector 
i = n x 7,  which  will  be  called  the  (unit)  cross-boundary  tangent  vector.  It  shall  be 
referred  to  as  if  it  were  provided  by  the  user.  We  will  calculate  a filling  surface  patch 
that  interpolates  the  user  supplied  boundary  curves  and  has  the  same  surface  normal  in 
a number  of  points.  This  will  leave  us  some  degrees  of  freedom,  which  we  will  use  to  fit 
the  curve  and  the  cross-boundary  tangent  vector  in  between  each  pair  of  interpolation 
points.  In  section  2 we  briefly  recall  the  basic  properties  of  uniform  Powell-Sabin  splines. 
Section  3 explains  how  we  can  benefit  from  these  properties  to  use  UPS-splines  for  the 
polygonal  hole  problem.  Section  4 explains  our  algorithm  in  detail.  Finally  we  remark 
that  on  the  pictures,  we  will  denote  2D  and  3D  entities  interchangebly;  therefore  most 
pictures  reflect  the  situation  only  schematically. 
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Fig.  1.  User  supplied  data. 

2 Uniform  Powell-Sabin  splines 

This  section  recalls  the  main  properties  of  Uniform  Powell-Sabin  splines.  For  details, 
we  refer  to  the  original  papers  [3,  5]. 

By  $2  (A*)  we  denote  the  linear  space  of  uniform  Powell-Sabin  splines  (in  the  sequel 
called  UPS-splines),  i.e.,  piecewise  quadratic  polynomials  on  a uniform  triangulation  A 
(which  means  that  all  triangles  are  equilateral  and  have  the  same  size)  of  a polygon  Cl, 
where  A*  is  a PS-refinement  of  A.  The  boundary  of  Q will  be  called  SCI,  whereas  the 
boundary  of  the  triangulation  will  be  referred  to  as  SA.  The  vertices  of  A are  denoted 
Vi,i  = 1 and  its  triangles  are  pt,i  — 1 These  splines  have  global  C1- 

continuity  on  A*.  Any  s(u,v)  has  a unique  B-spline  representation 

n 3 

s(u,v)  = ^^CijBJi(u,v),  (u,v)eO,  (2.1) 

i=i j=i 

where  the  locally  supported  basis  functions  form  a convex  partition  of  unity  and  cy  S R3 
are  the  control  points.  It  follows  that  s(u,v)  belongs  to  the  convex  hull  of  {cy  .. 
Furthermore,  one  can  prove  that  the  control  triangles,  being  defined  as  T^c^i, c^a, Cj,s), 
i — 1 ,...,n,  are  tangent  to  the  surface  at  s(Vi).  Due  to  the  local  support  of  Bf,  a 
change  to  cy  will  only  affect  s(u,  v)|mp  i.e.,  the  restriction  of  s(u,v)  to  the  molecule 
of  Vi,  being  the  set  of  triangles  p0  that  have  V)  as  a vertex.  This  indicates  that  we 
have  a useful  representation  for  C1-continuous  surfaces,  without  being  restricted  to  a 
rectangular  domain,  and  still  enjoying  the  interesting  features  of  the  classical  B-spline 
representation  for  tensor  product  splines. 

2.1  Subdivision 

In  [5]  we  present  a subdivision  scheme  for  UPS-splines.  Let  Ar  be  a uniform  refine- 
ment of  A,  obtained  by  midedge  subdivision.  For  a given  s(u,  v)  on  A,  the  representation 
(2.1)  on  Ar  can  be  calculated  using  convex  barycentric  combinations  of  the  control  points 
only.  First,  a new  control  triangle  along  each  edge  VfVj  is  calculated  as  illustrated  in 
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Fig.  2.  Subdivision  and  Bezier  points. 

Figure  2,  left,  for  the  bottom  edge  of  a triangle  p;(Vri,  V,  , 14)  G A: 

Cl  = 5(ci,2  + Ci,3) 

c2  = Icj,i  + i(ci,2  + Cj,2)  (2.2) 

C3  = Icj,l  + |(Ci,3  + Cj,3). 

Next,  the  control  triangles  at  the  original  vertices  are  rescaled:  for  example, 

Ci,l  = |ci,l  + g(Ci,2  + Ci,3) 

Cl,2  = |cI,2  + g(Ci,3  + Ci,i)  (2.3) 

Ci,3  = |ci,3  + |(Ci,l+Cii2). 

They  are  still  tangent  to  the  surface  at  their  barycenter,  but  their  area  is  only  a quarter 
that  of  the  former  control  triangles.  Therefore  they  connect  tighter  to  the  surface. 

2.2  The  piecewise  Bezier  representation 

Another  important  property  of  the  B-spline  representation  for  UPS-splines,  is  that 
the  piecewise  Bezier  representation  can  be  calculated  from  (2.1)  using  simple  convex 
barycentric  combinations  of  the  control  points.  In  particular,  focus  an  edge  F,V)  of  A 
(see  Figure  2,  right).  The  Bezier  points  of  the  edge  curve  can  be  found  from: 

s(Vi)  = Pi  = |(ci,i  + Ci)2  + cii3),  s(Vj)  = pj  = i(cj>i+cj,2  + cj|3),  (2.4) 

1 2 1 1 

Ui  = -(cii2  + Ci)3),  Uj  = -Cj,i  + -(cjj2  + Cj,3),  r,j  = -(Ui  + Uj).  (2.5) 

This  is  a piecewise  quadratic  Bezier  curve,  which  means  that  pi,  rjj  and  pj  are  surface 
points,  and  that  iq  - p;  and  pj  - Uj  are  tangent  to  the  surace  at  p;,  resp.  pj.  Assuming 
a (counterclockwise)  ordering  of  the  boundary  vertices  V)  € <5 A,  the  edge  curve  from 
s(Vj)  to  the  next  adjacent  point  s(Vj)  will  be  denoted  ei(u,  v). 

3 Application  to  the  polygonal  hole  problem 

Recall  that  our  goal  is  to  calculate  a UPS-spline  filling  a hole  in  a surface,  given  by  a 
set  of  bounding  curves  (denoted  p),  their  derivatives  7 and  the  cross-boundary  tangent 
vectors  5.  The  UPS-patch  will  fit  these  curves  approximately  along  its  boundary.  In  the 
first  place,  interpolation  of  the  given  data  at  the  vertices  Vi  G <5  A is  achieved.  This  leaves 
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Fig.  3.  Tangent  and  cross-boundary  tangent  vectors. 


some  degrees  of  freedom  allowing  to  fit  the  given  curves.  In  the  sequel  we  shall  denote 
the  user  supplied  data,  evaluated  at  V),  by  (pi,7i,<$i). 


3.1  Interpolating  UPS— splines  and  degrees  of  freedom 

In  order  to  obtain  interpolation  we  determine  a control  triangle  T*  in  the  tangent  plane 
spanned  by  pi  + ey  + i/S-,,  e,  v € R,  such  that  s(Vi)  = pi.  Curve  point  interpolation  is 
simply  expressed  by  (2.4).  Furthermore,  we  let  the  tangent  to  ej  at  Vi  be  parallel  to 


ui  - Pi  = -(ci,2  + Ci,3)  - -Ci,i  = ai% 


(3.1) 


where  a,  is  a scaling  factor.  Next,  we  need  the  cross-boundary  tangent  vector  of  s(u,  v) 
at  Vi  to  be  parallel  to  <5*.  Mapping  the  cross-boundary  vector  d in  the  domain  plane  (see 
Figure  2,  right)  onto  the  control  triangle  yields  a vector  parallel  with  Ci>2  — Cij3: 


Ci,2  ~ ci>3  = 2M, 

where  Pi  is  again  a scaling  factor. 


(3.2) 


Solving  (2.4),  (3.1)  and  (3.2)  to  Cy  in  terms  of  the  unknown  a,  and  /?*  (further  called 
the  a-  and  /3-factors)  yields 

{Ci,l  = Pi  - Qifi 

ci,2  — Pi  + ^7i  + Pi5. i (3.3) 

Ci,3  = Pi  + ^7i  - P&- 


These  equations  ensure  that  s(u,  v)  interpolates  the  given  data  at  Vi  e 6 A,  and  leaves 
us  two  degrees  of  freedom  per  vertex  (a*  and  pi).  These  scaling  factors  are  related  to 
the  size  of  the  control  triangle.  For  example,  subdivision  by  (2.3)  divides  a*  and  Pi  by  a 
factor  of  2. 


3.2  The  fitting  equations 


We  will  now  use  these  degrees  of  freedom  to  fit  the  user  supplied  data,  in  between  each 
pair  of  adjacent  interpolating  vertices  V.,  V)  € S A.  First,  the  a-factors  at  V)  and  Vj  are 
determined  by  trying  to  interpolate  the  curve  p at  the  edge  midpoint  V,j  = |(Vi  4-  Vj). 
From  Section  2.2,  the  interpolation  condition  reads  ry  = |(ui  + uj)  = py,  where  py 
is  the  given  curve  point.  Taking  (2.5)  and  (3.3)  into  account,  we  have 

aili  ~ <*jlj  = 4py  - 2(pi  + pj)  = qy.  . (3.4) 
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This  is  a system  of  3 equations  with  (at  most)  2 unknowns.  It  can  be  solved  in  the  least 
squares  sense. 

Next,  the  /3-factors  at  Vt  and  Vj  are  obtained  by  fitting  the  cross-boundary  tangent 
vector  at  Vjj.  First,  we  derive  a subdivision  rule  for  the  /3-factors  at  the  vertices  of  A 
from  (2.2)  and  (3.2): 

&A,j  = + M).  (3-5) 

where  5[  • is  the  cross-boundary  tangent  vector  to  s(u,  v)  at  Vij.  This  fiV- factor  belongs 
to  a finer  subdivision  level  then  /3,  and  /3 j,  so  we  have  to  scale  it  up  by  a factor  of  2.  The 
interpolation  condition  then  is 

Pijhj  = i(M  + PjSj).  (3.6) 

Note  that  Sij  has  been  used  instead  of  S't  J . This  is  again  an  overdetermined  system 
which  can  be  solved  in  the  least  squares  sense. 

4 The  algorithm 

We  will  restrict  the  figures  illustrating  the  algorithm  to  the  case  of  a triangular  hole, 
although  the  algorithm  is  immediately  applicable  to  cases  with  4,  5 and  6 boundary 
curves  as  well  (see  Section  4.4). 


The  idea  is  to  calculate,  during  a pre-iteration  step,  an  initial  solution  which  is  smooth, 
but  in  general  not  close  enough,  and  to  refine  this  approximation  iteratively  to  obtain 
a better  fit  to  the  given  curves  until  a certain  stopping  criterion  is  satisfied.  Finally, 
during  a post-iteration  step,  the  interior  control  triangles  are  calculated,  actually  filling 
the  hole.  Figure  4 illustrates  this:  imagine  a pre-iteration  step,  two  refinement  steps  and 
a post-iteration  step.  The  control  triangles  added  during  a particular  step  have  been 
shaded. 


4.1  An  initial  solution 

The  initial  solution  (Figure  4,  leftmost)  is  easily  obtained  by  solving  (3.4)  in  the  least 
squares  sense  for  each  edge  V)V).  If  we  assume  that  y,  ^ 7 ),  then 

: ..  a*  = ^ ((7.  ’ qy)  - (7j ' Qij)(7i ' 7j))  > 


(4.1) 
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ai  = ^ (~(7j  ' Qij)  + (7i ' qy)(7i ' 7j)) , (4.2) 

where  D = 1 — (7, -7j)2.  This  yields  two  o-factors  per  vertex:  one  for  each  boundary  edge 
being  incident  to  that  vertex.  Therefore,  2)  is  completely  determined.  The  /3-factors  can 
be  calculated  by  writing  (3.3)  for  both  edges  incident  with  the  vertex  and  eliminating 
c2,  respectively  Ci,  e.g.,  for  Figure  3,  right, 

/3i  = 0:2(72  • <5i),  P2  = -0:1(71  • ^2)-  (4.3) 

There  exist  pathological  cases  where  72  X <5i  or  71  X 62.  Our  algorithm  then  sets 
Pi  = Oi,  resp.  P2  — 02.  For  the  case  7 \ = 7 ),  (3.4)  has  no  solution  in  the  least-squares 
sense.  Assuming  that  si  is  a straight  line  from  s(Vi)  to  s(Vj),  the  o-factors  can  then  be 
determined  from  the  projection  onto  the  domain  plane,  where  the  size  of  the  so-called 
PS-triangles  (the  projections  of  the  control  triangles)  is  fixed.  The  reader  can  verify  that 
this  yields  o,;  = Oj  = 

4.2  The  iteration  step 

First  the  control  triangles  from  the  previous  steps  are  rescaled  by  subdivision.  This  is 
simply  done  by  scaling  down  the  a-  and  /3-factors:  a,;  <—  and  /3,  for  each 

Vi  € SA.  Next,  a new  control  triangle  is  created  in  between  any  two  adjacent  vertices  at 
the  coarser  level.  This  situation  is  illustrated  in  Figure  5,  left,  where  the  darker  triangles 
are  known.  We  are  looking  for  the  n-and  /3-factors  for  the  middle  control  polygon,  which 
is  tangent  to  the  surface  at  s(Vk),  V*  = |(Fi  + Vj).  Consider  the  a- factor  first.  In  order 
to  obtain  a better  fit,  we  try  to  interpolate  p at  Vi)k  = |(Vi  + Vk)  and  Vk.j  — 5 (Vk  + Vj). 
This  yields  a set  of  fitting  equations 

/ — aklk  = 

\ akik-atfj  = qkjj, 

where  a,  and  otj  are  known.  Thus,  ak  can  be  obtained  as  the  least-squares  solution  of 

(4-4):  . 

ttfc  = 2 (Tfc  • (arfi  - TiX  + Qkj  - Q!j7j)).  (4.5) 


The  Pk -factor  is  found  by  fitting  the  cross-boundary  vectors  at  Vi}k  and  Vk.j,  i.e.,  by 
solving  the  following  system  in  the  least-squares  sense: 


Pi,k$i,k  — \{Pi&i  + PkSk), 

Pk.jbkJ  = 2 X pj  bj  ) , 


(4-6) 


where  Pi  and  pj  are  known.  If  5ijk  = 5k  = as  is  always  the  case  for  a planar  curve, 
this  system  has  no  solution  in  the  least-squares  sense.  The  pk  factor  can  then  easily  be 
obtained  by  equation  (3.6),  i.e.,  by  subdivision  and  upscaling. 


4.3  The  interior  control  points 

Finally,  as  soon  as  the  user  supplied  edge  curves  have  been  approximated  well  enough, 
the  interior  control  points  at  the  eventual  refinement  level  have  to  be  calculated.  We  will 
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Fig.  5.  The  refinement  and  post-iteration  steps. 


Fig.  6.  The  hole  and  the  triangular  patches. 


discuss  three  possibilities  by  the  help  of  an  example;  Figure  6 shows  a hole  (left)  and 
two  filling  patches  (right). 

Copy  From  Initial.  The  interior  control  points  are  obtained  directly  from  the  initial 
solution  by  subdivision.  This  guarantees  that  the  interior  of  the  patch  is  smooth.  A 
disadvantage  is  that  the  inner  of  the  first  approximation  in  general  has  no  connec- 
tion with  the  shape  of  the  edge  curves.  This  can  cause  unwanted  artefacts  near  the 
boundary,  after  a few  iterations  (see  Figure  7,  left).  The  next  option  will  therefore 
take  edge  features  into  account. 

Averaging.  We  will  fill  the  hole  gradually  by  calculating  a ring  of  control  triangles 
during  each  pass,  going  from  the  edge  towards  the  inner  of  the  patch.  Figure  5,  right 
shows  an  example  where  each  ring  has  a different  shade  of  grey.  At  each  step,  a 
control  triangle  of  the  current  ring  is  obtained  by  averaging  six  surrounding  control 
triangles.  These  come  from  the  initial  solution,  or,  if  possible,  from  a previously 
calculated  ring.  Edge  features  are  now  smoothed  out  towards  the  inner  of  the  patch. 
However,  there  is  a main  disadvantage  to  this  approach,  if  averaging  is  applied  after 
the  last  iteration  step:  the  unwanted  artefacts  mentioned  before  are  now  repeated 
for  every  ring,  smoothed  out  towards  the  inner  of  the  surface,  as  shown  on  Figure 
7,  middle. 

Instant  Update.  A good  compromise  would  be  to  take  edge  features  into  account 
before  we  finish  iterating.  This  can  be  accomplished  by  subdividing  the  initial  solu- 
tion at  each  refinement  step,  but,  we  always  overwrite  its  edge  with  the  most  recent 
boundary  approximation.  The  results  of  this  strategy  are  depicted  in  Figure  7,  right. 

In  any  case  can  the  user  change  the  interior  control  triangles,  and  still  he  has  a C1- 
continuous  filling  patch,  fitting  the  specified  edge  curves  with  demanded  precision. 
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Fig.  7.  Copy  from  initial  solution  and  averaging  (4  iterations);  instant  update  (3  iter- 
ations). 


Fig.  8.  Cases  with  4,  5 and  6 boundary  curves. 


4.4  A note  on  the  number  of  edges 

The  algorithm  sketched  in  Section  4 is  immediately  applicable  to  problems  with  4,  5 
and  6 boundary  curves  as  well.  Figure  8 shows  the  configuration  of  the  initial  solution 
for  each  of  these  cases.  If  we  are  working  with  5 edges,  there  are  2 edges  having  a 
control  triangle  at  its  midpoint  (shaded  darker).  This  requires  a tiny  modification  to  the 
calculation  of  the  initial  solution  for  those  edges.  The  a-factors  are  obtained  by  solving 
(4.4)  to  the  unknown  , ay  and  a*.  The  /3-factors  of  the  outer  control  poygons  are 
obtained  as  usual;  for  the  middle  polygon  one  can  apply  (3.6).  Also,  for  the  cases  of  5 
and  6 boundary  curves,  an  interior  control  triangle  (unshaded)  has  to  be  calculated  for 
the  initial  solution.  This  can  be  done  by  averaging  the  six  surrounding  control  polygons. 
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